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Abstract 

We compute the normal forms for the Hamiltonian leading to the epicyclic 
approximations of the (perturbed) Kepler problem in the plane. The Hamil- 
tonian setting corresponds to the dynamics in the Hill synodic system where, 
by means of the tidal expansion of the potential, the equations of motion take 
the form of perturbed harmonic oscillators in a rotating frame. In the unper- 
turbed, purely Keplerian case, the post-epicyclic solutions produced with the 
normal form coincide with those obtained from the expansion of the solution 
of the Kepler equation. In all cases where the perturbed problem can be cast 
in autonomous form, the solution is easily obtained as a perturbation series. 
The generalization to the spatial problem and/or the non- autonomous case 
is straightforward. 

Keywords: Celestial mechanics. Planets and satellites: dynamical 
evolution and stability. Methods: analytical. 



1. Introduction 

The epicyclic approximation of the Kepler problem consists in the super- 
position of the motion of a guiding center along a circle and a motion along 
an ellipse of axis ratio 1:2. The radius of the main circle coincides with the 
semi-major axis a of the exact Keplerian orbit. The semi-minor axis of the 
ellipse is a factor e smaller, with e the eccentricity of the orbit. 

The power of the epicyclic formulation lays in being explicitly expressed 
in terms of the mean anomaly, therefore of time. This property is preserved 
by the 'post-epicyclic' approximation, namely the superposition of higher 
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harmonics of the mean motion providing a series of periodic deformations of 
the fundamental ellipse. 

The epicyclic and the post-epicyclic approximation can be considered as 
a perturbation expansion in which the circular orbit is the zero-order solu- 
tion, the epicycle is a first-order term and the post-epicycles are higher-order 
corrections. The eccentricity is the natural small perturbation parameter. 
A standard way to get the post-epicyclic approximation of the Keplerian 
orbits is through the series expansion of the Kepler equation as given by 
Laplace. There is an alternative approach to get solutions in the form of 
series by casting the problem in the perturbative setting of the Hill (or syn- 
odic) system: the dynamical problem is equivalent to that of an oscillator 
in a uniformly rotating frame. The linear harmonic approximation gives 
the standard epicyclic solution, nonlinear coupling terms give rise to higher 
harmonics. 

There has recently been a revival of interest in the post-epicyclic ap- 
proximations (the LISA experiment (Durandhar et al., 2005; Sweetser, 2005; 
Pucacco et al., 2010), coorbital motion of Saturn satellites (Gurfil and Kas- 
din, 2003), the J2 problem (Kasdin et al., 2005), etc.) due to the fact that 
the Hill frame provides a quite accurate description of perturbations of low 
eccentricity orbits, just by solving a simple quasi-linear mechanical system. 
It can be shown that, up to second order, the post-epicyclic approximation 
and the perturbed solution in the Hill system coincide (Rajesh Nayak et 
al., 2006). A natural question is up to which order, if any, this equivalence 
persist. In this regard, we recall here a classic result of celestial mechanics: 
the unperturbed Keplerian expansion is correctly reconstructed at arbitrary 
order by solving the equations of motion in the Hill system (Brouwer and 
Clemence, 1961). Although the same technique could be employed also to 
treat the perturbed case, we prefer to take a different approach. Aim of this 
paper is to investigate the perturbed Keplerian problem by constructing a 
Hamiltonian normal form and to get from it perturbation series solutions: we 
verify that the expansions for the unperturbed Keplerian orbits are correctly 
reconstructed. Then we proceed to solve samples of perturbed Keplerian 
problems. 

The plan of the paper is as follows: in Section 2 we recall the relation 
between the series solutions of the Kepler equation, the epicyclic formulas 
and the construction of the Hill system for the equations of motion; in Section 
3 we introduce the Hamiltonian normalization setting that we use to cope 
with the solution of the perturbed system in the synodic frame; in Section 4 
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we construct the normal form and compare its prediction with those obtained 
from the Kepler equation; in Section 5 we apply this formalism to the study 
of some examples of perturbed Keplerian problems; in Section 6 we give our 
conclusions. 



2. Epicyclic approximations of the Kepler problem 



We consider a Keplerian orbit of semi-major axis a and eccentricity e on 
the plane X, Y in the field generated by the mass M. A simple way to get the 
epicyclic approximation is by inserting the series expansion of the solution 
of the Kepler equation 

u — e sinw = nt, (1) 



into the parametric form 



X 
Y 



a[cosu — e 



e A sin u 



(2) 
(3) 



of the orbit. In these expressions u is the eccentric anomaly, the mean motion 
is 



n 



GM 



(4) 



and the product nt is referred to as the mean anomaly also denoted with I 
in the following. The series expansion in the eccentricity of the solution of 
the Kepler equation ([I]) is (Boccaletti and Pucacco, 1999) 



u = nt + e sin nt + -e sin 2nt + e 



sin 3nt 



sin rat 1 + 



(5) 



Inserting this into ^ and ^ and expanding in series of the eccentricity 
gives two Fourier series that can be referred to as epicyclic expansions of the 
Keplerian orbit: 



X 
Y 



cos nt — — (3 — cos 2nt) + 



sin nt H — sin 2nt + 
2 



(6) 
(7) 



In FigJTJ which is inspired by a picture appearing in Sweetser (2005), we can 
see the role of the different terms in the expansions f3|-[7|): the uniform motion 
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on the great circle is given by the trigonometric terms in nt of amplitude a; 
the motion on the ellipse is produced by the terms in 2nt of amplitude cte/2 
displaced along X by the constant amount — 3ae/2. The superposition of 
these two motions provides the first-order approximation of the Keplerian 
orbit. 




Figure 1: Epicycle representation of a Keplerian orbit. 



For our needs it is more convenient to express these expansions in the 
frame rotating with the mean motion (synodic frame) which turns out to 
be ideally suited to this perturbative approach. We make the change of 
coordinates 



so that 



x 



ae 



ae 



x = 

y = 



— cos nt — 



X cos nt + Y sin nt - 
Y cos nt — X sin nt, 



[1 — cos2nt) + -e 2 (cos3nt — cos nt) + ... 
8 



2 sin nt 



- sin 2nt H — e 
4 8 



sin 3nt — 3 sin nt 



(8) 
(9) 



,(10) 
(11) 
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The first-order terms in these expansions provide the properly called epicyclic 
solution: an ellipse with the semi-minor axis of length ae along the (radial) 
x-axis and the semi-major axis of length 2ae along the (negative, tangential) 
2/-axis. The higher-order terms are higher harmonic periodic corrections. In 
Figj2] we see the comparison between the exact Keplerian solution and the 
first two terms in the expansion: an unusually high value of the eccentricity 
is chosen for the plots (e = 0.5) in order to magnify the discrepancy between 
the true orbit and its approximations. 




20 40 60 80 100 



Figure 2: Comparison between the exact Keplerian solution (dots) of an orbit with a = 
1, e = 0.5 and the first two terms in the post-epicyclic expansion: the continuous curves 
represent the epicyclic (first-order) terms; the dashed curves include the contribution of 
the second-order terms. 

An alternative route to the epicyclic approximation is that based on the 
solutions of the perturbed equations of motion in the Hill frame which is 
convenient when we have to analyze in the simplest way the geometry of low 
eccentricity orbits and their most relevant perturbations. The Hill reference 
system (see, Boccaletti and Pucacco, 1999 and Bik et al. 2007), also called 
Clohessy- Wiltshire system (Clohessy and Wiltshire, 1960) allows us to reduce 
the equations of motion to those of a set of coupled oscillators. The Hill 
equations properly said consist of the linear part of this system which is 
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readily solved. Non-linear and non-autonomous terms can be added and 
solved perturbatively. 

The synodic reference system given by transformations ( 8p ) is based 
on a circular orbit of radius a around the central mass and has its origin 
rotating at angular velocity n with the reference orbit. In the applications 
to a perturbed Kepler problem, it is useful to relax this assumption and the 
angular velocity can be arbitrary, say Q. Reverting for the time being to 
the three dimensional description for a clearer geometric view, we have the 
axes x, y, z such that x is directed radially opposed to the center, y is in the 
direction tangent to the motion and z is perpendicular to their plane. The 
new equations of motion are introduced according to the following steps: 

1. Transform from the inertial barycentric system to the non-inertial 
synodic system: the equations of motion take the form 



x = -V$ - 2Q A x - Q A (SI A (a: + a)) . (12) 

Here the dot denotes the derivative with respect to time, $ is the total 
gravitational potential and x — (x,y, z), ft — (0, 0, ft), a — (a, 0, 0). 
2. The potential is approximated according to the expansion 



$ = $1 



x=0 



+ 



9$ 



dxy 



x\ + 



1 d 2 $ 



x=0 



2 dx^dx\ 



x^x\ + 0(|je| 



(13) 



x=0 



where the sum is implicit over the repeated indexes A,/i = 1,2,3. The 
"far tide" (or Hill) approximation (Boccaletti and Pucacco, 1999) consists 



of retaining only the quadratic terms in (13) and therefore only first order 



terms into the expansion of the gradient of the potential. 

Taking as reference example the unperturbed Kepler potential we have 



$ = $ 



GM 



K 



\/(x + a) 2 + y 2 + z 2 
so that the Keplerian far tide approximation gives 



<9$ 



K 



dx 



GM 2 <9$ K 

— 7T- — n a. 



x=0 



dy 



x=0 



dz 



(14) 



(15) 



x=0 



and 



dx 2 



2GM 



-2n 2 



x=0 



d 2 $ K 
dy 2 



x=Q 



d 2 $ K 
dz 2 



GM 



n 



(16) 



£C=0 
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all others zero. The next order of approximation (octupole) is: 



<9 3 $ K 




6GM 


Qn 2 


dx 3 


x=0 


a 4 


a 


<9 3 $ K 






3GM 3n 2 


dxdy 2 


x=0 


dxdz 2 


x =o « 4 a 



(17) 
(18) 



all others zero. By using these in (12), with Q = n, we get 



X 




2ny — 


3n 2 x + 






2nx — 


3n 2 


y 


+ 


— xy 






a 




+ 




3n 2 


z 


n 2 z — 


xz - 








a 



3n 2 

1 

2a 



(19) 



where in the vector f(x,e) we include the higher-order terms produced by 



the expansion (13) and, in general, autonomous perturbations. We assume 
that any other perturbing force can be derived by an analogous series expan- 
sion of the corresponding potential and added to this system. Each compo- 
nent, coherently with a perturbative approach, will be a power series in the 
coordinates with suitably 'small' coefficients e = (e x ,e y ,e z ). 

In these equations the second order approximation of the perturbed Ke- 
plerian motion is explicitly given. Coherently with the above treatment, we 
will limit the analysis to the motion on the plane and will drop the third 
equation for the z oscillation. This does not affect too much the generality 
of the approach and can be easily included if desired. 

The perturbation vector, in its simplest form, consists of constant and 
linear (in x) terms of the form Eq + E\X + . . . : no terms in y are allowed 
if we want preserve the autonomous nature of the perturbation. This gives 
explicit small changes already at zero-order. These and possibly higher- 
order perturbing terms propagate in the full solutions through the non-linear 
coupling with the other components of the field expansion. The presence of 
constant terms in the perturbation suggests to generalize the Hill equations 
by not specifying in advance the angular velocity and leaving it as the free 
parameter Q that will typically be a small perturbation of the Keplerian value 
Rather we leave a as the fixed radius of the reference circular orbit. In 
this way, the relevant system at first order can be written in the form 
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x - 2tty - tt 2 p x + 0(\x\ 2 ) = 



y + 2Qx + 0(\x\ 2 ) = 0, (20) 

where 

Q 2 p ~ (3 + 2e)fi 2 (21) 

includes the first order contribution of the perturbing field and the small 
parameter e is determined by the nature of this field. In this work we will 



consider the simple class of the form $ pcr t ~ 1/r 3 , but eq.(21) holds for any 



third-order expansion complying with the above stated conditions. Since in 



general the presence of quadratic or higher-order terms in the system (20) 
makes it non-integrable, the only viable route to get some solutions is to 
head to a perturbative approach. The solution of the system can therefore 
be written in the form of a series expansion 

x(t) = x w + e ■ x {2) + . . . , (22) 

where the value of the upper index denotes the order. 

The first-order general solution is obtained by solving the linear part of 



(20). The characteristic equation of this linear system is 



A 4 + A 2 (4ft 2 -^) =0, (23) 
so that one of its eingenvalues vanishes and the other provides the frequency 



io = J AVI 2 - VI 2 ~ n(l - e). (24) 



The zero eigenvalue is associated to the neutral stability of one of the normal 
modes therefore producing a linear drift. Choosing initial conditions in order 
to remove constant and secular terms due to this drift, the solution of the 
linear system is: 

= -Acos(ut + (j)), y {1) = 2-Asin(c^ + 0). (25) 

In the purely Keplerian case, 



to 



Q — n, Q p = \/3n, A = ea, = 0, (26) 
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they coincide with the first terms in ( TofTT ) and correspond to harmonic 
oscillations around the 'reference' x = circular orbit. In this case, the 
perturbative setting is provided by the observation that the radial amplitude 
is constrained by the apsides and therefore we can rescale the coordinates 
according to 



x 



x 



ae 



y 



y 



ae 



(27) 



Additional perturbations add gradually smaller terms according to the size 
of coupling constants like those above. Higher order terms in the solution 
can be obtained by plugging the series (22) in the system (20), equating 



terms of the same order in e and exploiting solution (25). On the other 



hand, a more general approach to the solution of perturbation problems can 
be obtained with the more effective method of constructing a Hamiltonian 
normal form, as shown in the next section. From the analytical point of view 
it is worth mentioning that even the simple system (19) with just f(x, e) = 
(with x(t) = x^ + ex^ as a second-order approximate solution) is not 
integrable, unlike the original Kepler system which is clearly integrable. Is 
a common lesson of modern analytical mechanics that truncated solutions 
of non-integrable approximations can be more useful than formally exact 
solutions of their integrable counterparts. 



3. The Hamiltonian setting 

The equations of motion in the Hill system can be obtained by introduc- 
ing a suitable Hamiltonian. Applying the approach usually adopted with 
Hamiltonians admitting linear terms in the momenta (Hietarinta, 1987; Pu- 
cacco and Rosquist, 2005) it can easily be proven that the planar [z = 0) set 



of system (19) is generated by the Hamiltonian 



U{p, x) = l -{ Vx + nyf + \{p y - Slxf - \ny + <&(*>, (28) 

where $^ 3 ^ is the expansion of the potential (Keplerian and, possibly, per- 
turbations) starting from the cubic terms. The canonical equation given 



by this Hamiltonian coincide with the first two equations of system (19) 
if we retain in the expansion of $k only the cubic terms and assume that 
f(x,e) = -V§ pert (x,£). 
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The most effective way of investigating the solution of the canonical equa- 



tion generated by (28) is by constructing a normal form by a suitable canon- 
ical transformation. We look for a new Hamiltonian in the new canonical 
variables P, X , given by 

oo 

/C(P,X) = ^/C n (P,X), (29) 

n=0 

with the prescription that 

{G,JC} = (30) 

where Q is a suitable low degree polynomial in the canonical variables. When 
Ho, the quadratic part of the Hamiltonian under investigation, is semi-simple, 
the usual choice is simply Q = H . In the present case this is not possible 
and we have to make a different choice. In any case, in order to set up the 
procedure based on the Lie transform method, %q must be diagonalized (in 
these and subsequent formulas we adopt the convention of labeling the first 
term in the expansion with the index zero: in general, the 'zero order' terms 
are quadratic homogeneous polynomials and terms of order n are polynomials 
of degree n + 2). 

With the canonical transformation (p x ,p y ,x,y) — > (P£,Pr],£,,v) defined 

by 



we get 



where 



(31) 



y = jwte-v), (32) 

p x + tty = y/upt, (33) 
p y -ttx = -^( ^^ Pv-Wn, (34) 



-Ho = l -u{p\ + £ 2 ) + \k V \ (35) 



4^ - ui + ni' [6b) 
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This quadratic part is of semi-simple+nilpotent type: this case has already 
been studied (Broer at al., 1993; Palacian and Yanguas, 2000) and it appears 
that a natural choice to work out the normalization is to extend only the 
partial energy 

G = \{$+e\ (37) 

which is an obvious integral of motion of Hq, to the whole normal form /C. 

In the process of normalization based on the Lie transform, new coordi- 
nates P, X result from the canonical transformation 

(pe,t) = T x (P,X). (38) 

Considering a generating function %, the Lie transform operator T x is defined 
by (Boccaletti and Pucacco, 1999) 



T x ^M fc ( 39 ) 



fc=0 

where 

k . 

M = l, M k = Y d 3 -L Xj M h . j . (40) 
i=i 

The functions Xk are the coefficients in the expansion of the generating func- 
tion of the canonical transformation and the linear differential operator L x 
is defined through the Poisson bracket, L x {-) = {x, •}■ 

The terms in the new Hamiltonian are determined through the recursive 
set of linear partial differential equations (Giorgilli, 2002) 

/Co = 'Ho, 

/Ci = Ux + MiQ = H! + L X1 Q, 

jc 2 = u 2 + m 1 u 1 + M 2 g = u 2 + l x1 u x + \l 2 xi u g + L X2 g, (41) 



J 3 
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where £?.(•) = {Xi {x> '}}• 'Solving' the equation at the k-th step consists of 
a twofold task: to find JCk and Xk- Once a given 'target' order N = fc max has 
been reached the process is stopped and two truncated series 



N 



= J2xk(P,X), (42) 

fe=i 
v 

= J2)C k (P,X), (43) 



fc=0 



are constructed. The new Hamiltonian /C( fcmax ) is integrable because, by con- 
struction, admits an additional integral of motion. If we are able to solve 
canonical equations, we get the flow (P(t), X(t)) and we can use the generat- 
ing function ^( fcmax ) to construct the transform (38) and the transformation 



rules (31-34) to obtain the flow in terms of the initial 'physical' variables 
(p{t), x(t)). Also these expressions will be series truncated at the same order 
as the generating function. 

In the following section we apply this theory to the Hamiltonian of the 
Keplerian problem in the Hill frame by explicitly constructing the series 
x(t) = + ex^ + ... showing that they coincide with that obtained by 
expanding the series solution of the Kepler equation. Later on we outline the 
general approach to cope with autonomous perturbations. 

4. Normalization and construction of the generating function for 
the unperturbed Kepler problem 



We have explicitly constructed the two truncated series (42) and (43) 
in the Keplerian case up to A; max = 4, namely up to and including terms 
of degree 6 in the canonical variables. We start by expanding the Kepler 
potential up to order 6, giving the functions 
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n ( 3 3xt/ 2 



71 / 4 , o 2 2 



H 3 = (x 5 -5xV + , (46) 



a 3 



n ( a 15x 4 w 2 45x 2 w 4 5y 6 . 
^4 = "7 \-^ + ^ L -^ L + ^r • (47) 



™ r_3..2 , 15x ?/ 4 

3 

a 1 V ~ ' 2 8 '16 

These expressions are written with a small abuse of notation with respect to 
the previous section, where the higher-order terms of the perturbed Hamil- 
tonian are assumed to be already prepared for the normalization and should 



be expressed in terms of the new variables introduced by (31-34). However, 
it is more simple and natural to use the old physical coordinates. Moreover, 
the Hamiltonian is rescaled dividing by the mean motion so that the new 
time associated to the transformed Hamiltonian is the mean anomaly I = nt. 



Using (36) and the first two equalities of (26), the unperturbed diagonal part 
is therefore 

U Q = \{pl + e)-\pl (48) 

A first notable outcome of the normalization procedure is that, up to N = 
fcmax — 4, all the perturbing terms in the Hamiltonian can be removed and 
therefore the normal form up to this order is simply 

^/Q, = ^(Pl + X 2 )-^, (49) 

with JCj = for 1 < j < 4. We may conjecture that this is true for the Kep- 
lerian expansion at all orders but, anyway, we can imagine that this property 
is in general lost as soon as one introduces an arbitrary extra perturbation. 

The extremely simple structure of the normal form is counterbalanced by 
the generating function which, in view of the its role of eliminating a huge 
number of terms at each step, is quite cumbersome. We limit ourselves to 
give the explicit expressions up to order 4 which will be enough to establish 
our result. We get 



13 



an X i = ~X 2 Y + ^ + X 2 P X - 6Y 2 P X + 5YP* 

+ 3XFP y - 9 -XP x P Y - ^ + ^Py, (50) 



a 2 n 2 x 2 = — 6XY 3 - X 3 P X + 15XF 2 P X - Y1XYP 2 X 

, 37 , 7F 3 P y 7 9 

+3XP 3 - ^X 2 YPy + — ^ + |^ 2 PxPy 

--F 2 ? x P y + -YPlPy - -P\Py 

4 12 x 12 x 



and 



1 1 X yp 2 - -xp x p y - ^2 - 1, 

8 y 16 y 6 32 



aVxs = ^ — + 



101X 4 y 200X 2 F 3 268F 5 7X 4 P X 



36 9 45 



6 3 18 x 9 x 



31 „ 283 725FP4 lOlPf 

X 2 Pl + F 2 P|- - + - 

3 x 6 x 36 30 

245 o 58 o 65 o 293 , 

+ — X 3 YP Y - -XY 3 P Y - ^^ 3 ^xPy + —XY 2 P X P Y 

U5 XYP 2 P Y + —XP X P Y - —X 2 YP 2 + ^F 3 P 2 



4 16 108 54 

+ ^X 2 P X P 2 - -F 2 P X P 2 + —YP 2 P 2 - —P 3 X P 2 
12 8 y 27 x y 144 x Y 

48 y 24 y 108 192 y 

We recall that, at each step of normalization, we get new canonical vari- 
ables on which the normal form depends. However, in view of the result 
(49), at each order the solution of the canonical equations given by the new 
Hamiltonian (49) 
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dX 

~dl 

dY 

~dl 



Px, 
3 

~4 J 



dP 



x 



-X: 



dPs 



Y 



Y; 



0. 



(51) 
(52) 



is the same. With a natural choice for the arbitrary constants we can elimi- 
nate secular terms and adjust amplitudes so to have the normal mode 



X nm (£) = A x cosi, Y nm (£) = (53) 

where Ax is an arbitrary amplitude at our disposal for any needed adjust- 
ment. To perform the comparison with the exact epicyclic expansion, we can 



now apply the rule given by the general transform (38) to the solution on 
the normal mode to find the series 



iO-) 



j(2) 



£ (3) + ... 



Pi 



(1) 



Pi 



(2) 



Pi 



(3) 



(54) 
(55) 



Transformations (31-34) will finally provide the series in the initial 'physical' 



variables. The transformations are expressed by the same relations as in the 



general procedure sketched in (41) and to second order, in addition to the 
obvious identities 



are explicitly given by 



Pi 



(i) 



Pi 



r(2) 

(2) 



\Pxtuni Pymxl)) 



L xi X = {Xi, x }, 
L Xl P = { X i,P} 



(56) 
(57) 



(58) 
(59) 



and 



Pi 



t(3) 
(3) 



L X2 X + 2 L *i X 



--{ X 2,X} + -{xi,{Xi,X}}, (60) 
{X2,P} + lixx,{XuP}}- (61) 
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Finally, by using transformations ( 31p2 ), we obtain the predictions in terms 
of the original variables. 

To reconstruct in the correct way the whole series we have to properly 



set arbitrary constants. With the choice made in eqs.(53) we have only to 
adjust the amplitude of the normal mode. If we write it in the form of a 
series in the eccentricity 



A 



x 



(62) 



we can choose the arbitrary coefficients aj to satisfy the initial conditions. 
The Kepler orbit we are looking for is such that the periastron is at t — 0. 
Therefore a suitable initial condition is x(0) = — ae,y(0) = 0: at each order 
these constraints determine the value of the ctj, j > 1 and, as a consequence, 
the terms in the expansions with the correct coefficients. With the choice 



a i 



-a, we start with the first-order terms 



x 



(i) 



— aecosj 
2ae sin£ 



(63) 
(64) 



as in (ToJTT). The constraint at periastron impose the choice a 2 = at second 
order and a% = —a/8 at third order, so that we get 



x {2) = ae 2 - (cos2£- 1) 



/ 2 ) = ae 2 -sin2£ 



,1 
4 



(65) 
(66) 



and 



x 



(3) 



ae — (cos3£ — cos^ 



,(3) 



ae 



24 



sin 3£ 



sin- 



(67) 
(68) 



Comparing with ( TofTT ) we see that the predictions offered by the normal 
form approach are exact up to order e 3 . Going even farther, we get, at fourth 
order, 04 = and 
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^ = ae4 (l^ COS4 ^^ COS2£ "6^)' (69) 
y (4) = ae 4 ^sin4£-^sin2^, (70) 

which agrees with results obtained by patiently solving the equations of mo- 
tion (see e.g. Brouwer and Clemence, 1961). 

5. Normalization of the perturbed Kepler problem 

In the previous section we have validated the normal form approach by 
exploiting a high order (N=4) normalization of the epicyclic approximation of 
the Kepler problem. We can then rely on a powerful method of retrieving in- 
formation on the dynamics of the whole class of systems given by autonomous 
perturbations of the Kepler problem. 



Referring to the perturbation parameter introduced in (21) and (24), 
two important examples are: the 'Schwarzschild' correction mimicking the 
relativistic precession given by 

_3r s 2GM 
£S = 2^ rs = ^ (71) 
and a gravitational quadrupole correction given by 

3 , /iT 2 



^ S 2*UJ (72) 

where R is the radius of the central body. In both cases the parameters r s 
and J2 x R 2 are 'small' in these units (Pucacco et al., 2010). 

Another potential of great relevance is the so called isochrone (Boccaletti 
and Pucacco, 1999) that in this framework can be approximated by 

ei = --, (73) 
a 

where b is a length-scale assumed small with respect to a. 

It happens that, in general, the normal form is more complex than in the 
case of the Kepler problem. In particular, it is no more possible to eliminate 
all terms with the canonical transformation. We discuss below the origin and 
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the implications of this phenomenon. In our derivation we stop at the first 
steps of normalization procedure to avoid cumbersome formulas. However, 
this is still a very useful starting point for applications; in fact, in the spirit 
of perturbation theory and in view of the smallness of all coupling constants 
introduced above, terms in the solutions coming from these perturbations 
can be linearly superposed to terms giving the 'unperturbed' solution of the 
previous section. 

The order at which both series have to be truncated depends on the 
required precision and on the relative size of the perturbation terms. For 
example, in the case of the orbits of LISA (a = 1UA, e = 0.01), if one is 
interested in investigating the Schwarzschild precession due to the Solar field 
(es — 10 -8 ), in principle the Kepler expansion should be truncated at order 
e 4 : however, in practice, in view of the differential motion on the members 
of the 'constellation' that gives the measurable effect, one can truncate at 
order e 3 (Pucacco et al., 2010). In any case, we see that the greatest part of 
the 'job' is done by the Kepler expansion and one step of normalization of 
the additional perturbation is enough for our purposes. 

Therefore, we limit ourselves to writing the explicit form of the first non- 
vanishing terms of the normal form and of the generating function. They 
are: 



JC 2 — n 



^j{Px + x 2 ) 2 + 3or 2 (p| + x 2 ) + ior 4 



(74) 



and 



-X, = -iX 2 Y + ^f- + ^X'P X - WY 2 P X + 5YP$ 
Comparing this last expression with the corresponding one in the Keplerian 



case (50) we can see the slightly different approach in the two perturbation 
settings: here the perturbation parameter appears explicitly in the generating 
function (and in the normal form) whereas in the Keplerian case the role of 
the perturbation parameter is played by e and is hidden into the coordinates. 

Afterwards, proceeding as in the previous section, we obtain the second- 
order solution to the equations of motion representing the effect of the field 
perturbation: 
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X 



to 



7 13 

aee | coswit H cos 2u\t 

6 6 



,to 



2ae 



[1 + e) sin uj\t + — e sin 2uj\t 



(75) 
(76) 



This solution is one order (in e and e) beyond those found in the descrip- 
tion of the precessional effects to which is sensitive the LISA constellation 
(Pucacco et al., 2010); however, in that work, the first order (in e) terms 
have a more complex structure because that solution was computed in a Hill 
frame rotating at the fixed value n of the unperturbed Keplerian orbit. The 
additional constant part in x(t) appearing in Pucacco et al. (2010), —ea/3, 
accounts for the change of radius of the reference orbit due to the slight 
change in the force. Here this term disappears because the different angular 
velocity of the Hill frame gives a centrifugal force which exactly balances the 
total active force due to the Keplerian + perturbing field. A constant term 
appears only at second order. In Figj3]we show an example of the use of the 
two approximations. 




/(years) 



Figure 3: Approximations of the Schwarzschild precession with the post-epicyclic formulas 
in the Hill system: the x coordinate is on the left, the y coordinate is on the right. The 
continuous lines represent the residuals from the exact numerical evaluation obtained with 



the second-order expressions ( 75 - 76 1 whereas the dashed lines are obtained with the first 
order approach of Pucacco et al. (2010). The orbital parameters are a = 1UA, 



0.01. 



Another important difference in these solutions is the presence of the new 



frequency U\. This occurs precisely because the second order term (74) in the 
normal form is non-vanishing. If we write the generalization of the normal 



mode (53) in the form 
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X nm {t) = V2Jcoswit, Y nm {t) = 0, (77) 

where J is an action-like variable, its frequency can be computed from the 
Hamiltonian as 

- - % < 78 > 



where the derivative is evaluated on the normal mode (77). It is easy to 
check that 

(79) 

In practice, this result has negligible consequence on the predicted amount 
of these small effects: the advance of periastron is still approximated with 
high accuracy by the formula 

Aip = 2ir n ~ U ~ 2tie [rad/orbit] (80) 
to 

because the extra correction is of second order in the perturbation parameter. 
However, the result is important from a formal point of view, since it accounts 
for the dependence of frequency on the amplitude of oscillation. This remark 
shows us the relevance of the non-vanishing second (and higher) order terms 
in the normal form: only in this way this general phenomenon may happen. 
In the degenerate case of the exact Kepler problem (and in the trivial case of 
the harmonic oscillator) both frequencies Q and u are constant and coincide 
with the mean motion, so that = and all orbits are closed. In this 
perturbative setting this happens only if the normal form coincides with the 
quadratic term at every order. 



6. Conclusions 

The results presented in this work provide a comprehensive setting to un- 
derstand some stimulating but unsettled results obtained in previous studies 
(Gurfil and Kasdin, 2003; Kasdin et al., 2005; Durandhar et a!., 2005; Rajesh 
Nayak et al., 2006; Sweetser, 2005; Pucacco et al., 2010). 

We have here tried to address the natural question about the limits of 
validity of those results. We adopted the approach usually followed in exploit- 
ing approximate invariants constructed in a high order perturbation theory 
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(Giorgilli and Locatelli, 2006) and we have checked, in the specific problem 
of the perturbed Keplerian motion, how the power series representing the 
approximate solution can be used to imitate the exact solution. Overall, the 
approach based on the normal form extends in many ways the results of works 
based on the Hamilton- Jacobi equation and perturbations thereof (Kasdin 
et al., 2005) and on the simple direct analysis of the non-linear equations of 
motion (Rajesh Nayak et al., 2006; Pucacco et al., 2010). 

In the framework based on the Hill reference system, we can state that a 
perturbative approach makes sense if the expansion of the total potential is 
performed up to a degree such that the series representing the post epicyclic 
approximation of the unperturbed orbits and the amount of the perturbation 
give errors of comparable size. These results also provide answers to the 
questions raised by Sweetser (2005) concerning the relevance of the post- 
epicyclic series in satellite problems. In the study of satellite formations on 
orbits with small eccentricity this approximation is extremely good (Pucacco 
et al., 2010) and can be also profitably used in the coorbital problem (Gurfil 
and Kasdin, 2003) and in the so called 'tethered' system (Sidorenko and 
Celletti, 2010). 

The generalization to the 3-dimensional case can be performed with lim- 
ited effort. At the linear level, the oscillation in the z direction is decoupled 
and proceeds with its own proper frequency. In the presence of a pertur- 
bation, the non linear coupling is managed through the normal form with 
standard methods (Giorgilli, 2002). Non-autonomous perturbations can be 
considered by extending the phase-space to include the time variable. In 
general, in 3 dimensions, the normalizing transformation no longer provides 
an integrable normal form. Usually, the improvement consists only in lower- 
ing by one the dimensionality of the problem. However, the investigation of 
existence and stability of normal modes proceeds in much the same way as 
above. 
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